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Using grand canonical Monte Carlo simulations, we investigate the percolation be- 
havior of a square well fluid with an ultra-short range of attraction in three dimension 
(3D) and in confined geometry. The latter is defined through two parallel and struc- 
tureless walls (slit-pore). We focus on temperatures above the critical temperature 
of the (metastable) condensation transition of the 3D system. Investigating a broad 
range of systems sizes, we first determine the percolation thresholds, i. e., the critical 
packing fraction for percolation rj c . For the slit-pore systems, rj c is found to vary with 
the wall separation L z in a continuous but non-monotonic way, r] c (L z — > oo) = rf c . 
We also report results for critical exponents of the percolation transition, specifically, 
the exponent v of the correlation length £ and the two fisher exponents r and a of 
the cluster-size distribution. These exponents are obtained from a finite-size analysis 
involving the cluster-size distribution and the radii of gyration distribution at the 
percolation threshold. Within the accuracy of our simulations, the values of the crit- 
ical exponents of our 3D system are comparable to those of 3D random percolation 
theory. For narrow slit-pores, the estimated exponents are found to be close to those 
obtained from the random percolation theory in two dimensions. 
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I. INTRODUCTION 



Percolation is a geometrical transition, in which interacting units such as the particles in 
a fluid, the spins on a lattice, or the nodes of a network spontaneously form system-spanning 
clusters termed the "percolated phase".- In contrast the units are distributed homogeneously 
or form isolated clusters of finite size in the non-percolated phase. 

In this publication we focus on continuous systems, where percolation was originally dis- 
cussed as a phenomenon in the context of flow through porous media. Later, the statistical 
theory of percolation has been used to understand the critical behavior of fluid (Fisher's 
droplet model) particularly the divergence of the correlation length and the critical ex- 
ponents characterizing the behavior close to the critical point. Indeed, it is well known 
that percolation obeys concepts of universality and scaling, similar to second-order ther- 
modynamic phase transitions. 4 This includes the universality class of the so-called random 
percolation, which is assumed to apply to continuous systems such as fluids. 

More recent (experimental and theoretical) research on continuous percolation often in- 
volves colloidal suspensions. One important topic in this area concerns the percolation 
of rod-like colloids, prominent examples being carbon nanotubes and other carbon-based 
particles^-. From an applicational point of view the underlying idea is that the perco- 
lated network leads to lightweight materials with strongly enhanced mechanical stability 
and electrical (and/or thermal) conductivity. Another main topic, which concerns partic- 
ularly complex colloidal mixtures^ 1 ^ and colloids with directional interact ions^ 1 - - —, is the 
intimate relation between percolation and the formation of a physical gel, which is a state 
in which particles are connected via bonds of limited lifetime. It is now well established 
that such colloidal gels, which are characterized through a very specific dynamic behavior, 
can form at extremely low packing fractions. Note, however, that gelation is a phenomenon 
which normally occurs for very strong coupling conditions, i. e., at temperatures far below 
those related to the vapor-liquid critical point (if the latter exists at all). 

In the present study we are interested in the characteristics of the percolation transition 
at moderate (supercritical) temperatures, focusing on a system of spherical colloids with 
attractive interactions. Specifically, we aim to determine the percolation threshold, the 
scaling behavior and the related critical exponents. We consider both, three-dimensional 
(3D) systems and systems in slit-pore geometries which are infinite in only two dimensions 
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(2D). Our investigations are based on grand-canonical Monte Carlo (GCMC) computer 
simulations combined with a finite-size scaling analysis. 

Our study is partially motivated by a series of recent Monte Carlo (MC) results by 
Nezbeda and coworkers.— ~— These authors investigated the percolation in the corresponding 
supercritical regimes of various 3D (bulk) continuous model systems, differing in the precise 
form of the interaction potential.— Using different geometric cluster definitions and spanning 
rules, they found perfect scaling behavior for all systems investigated. However, the critical 
exponents turned out to be strongly dependent on the specific interaction potential and the 
temperature. This obviously contradicts the expectation that the percolation transition in 
these systems is universal. Indeed, in ideal random percolation the critical exponents and 
scaling functions only depend on the dimensionality and the symmetries of the system.- 

Here we present results from a GCMC study involving colloids with an ultra-short range 
of attraction, modeled by a square-well potential with an attraction width of only four 
percent of the colloidal diameter. This is considerably smaller than the range of fifty percent 
studied in Ref. Ij]- Interestingly, our GCMC results for the 3D case do indicate universality 
of the exponents in the sense that they do not depend on the temperature in the range 
considered. The exponents do, however, depend on the spatial dimensionality, consistent 
with the predictions from random percolation theory. 

We note that our model is by no means artificial. Indeed, colloids with ultra-short ranged 
interactions can be experimentally realized by adding small, non-adsorbing polymers to the 
colloidal solution, yielding short-range, attractive depletion interactions of tunable ranged. 
Another, recently studied example (where the justification of the ultra-short ranged square- 
well model has actually been checked using scattering data) are silica nanoparticles with 
added lysozym©^. 

This paper is organized as follows. In Sec. [Til we introduce the model and discuss its bulk 
phase diagram. Numerical details of the GCMC simulations are described in Sec. IIHI In 
Sec . 11 VI we present our numerical results for the percolation thresholds and critical exponents 
of the bulk and slit-pore systems. Specifically, we consider the exponents of the correlation 
length, the cluster size distribution, and the radius of gyration. In Sec. IIVCI we finally 
present results for the confinement-induced shift of the percolation threshold, and we provide 
a simple theory based on scaling arguments to describe this shift. We close this paper with 
the conclusions in Sec. [V] 



II. MODEL 



In our colloidal model system, the particles interact via an attractive square-well (SW) 
potential defined as 



OO, Tij < O 



a < m < A • (!) 



A < r. 



In Eq. ([!]) = \rij\ = \ri — Tj\ is the particle distance, a is the hard-core diameter, e > 
is the strength of the attraction and A determines the range of this attraction, i. e., A > a. 
Here we choose a value A = 1.04a, corresponding to a system with an ultra-short range of 
attraction. 

The SW model with ultra-short ranged attraction has attracted growing attention, since 
it is known to be a simple but adequate model for colloidal suspensions. Indeed, short-ranged 
attractive colloids can be realized experimentally by adding small non-adsorbing polymers 
to the colloidal solution. This creates a depletion effect leading to an effective attraction be- 
tween the colloids. To screen the remaining (attractive) Van der Waals interactions between 
the colloids, one chooses a solvent of similar dielectric permittivity.- 1 ^ 1 ^ If the difference in 
the length scales between the polymers and the colloids is sufficiently large, the degrees of 
freedom associated with the polymers can be "integrated out" of the theoretical description. 
The remaining implicit effect is an effective short-ranged attraction between the colloids 
due to depletion interactions.— Regarding the actual value of A, we note that for values 
A < 1.25a the stable liquid phase disappears and a typical colloidal phase diagram involving 
only a fluid and a solid phase arises.— 

We consider our colloidal model in three dimensions (bulk) and in slit-pore geometries. 
The latter are realized by two plane-parallel structureless walls, which are parallel to the 
x-y-plane of the coordinate system and located at ±L z /2. The wall-particle interaction is 
modeled via 

{oo, d < a/2 
( 2 ) 
0, d>a/2 

with the particle-wall distance d. Thus, the slit-pore-particle interaction is given by 

^siit(n) = n wall (L z /2 + zi)+ u WSL n(L z /2 -Zi)) (3) 
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where Zi is the z-component of the center of mass position T{ of the ith particle. Such 
a purely repulsive colloid-wall potential can be realized experimentally by using surface 
structures like a polymer-coated brush, which is penetrable for the polymers but not for the 
colloids.^ 

Recent research has shown that short-ranged systems, such as the present SW model, ex- 
hibit thermodynamic properties which are insensitive to the specific shape of the interaction 
potential and approximately fulfill an extended law of corresponding states.— ~— This will 
be useful to derive a rough draft of the bulk phase diagram of our system. To this end we 
first consider a limit case of the SW model, the Baxter model, which is also referred to as 
the adhesive spheres (AHS) model. Within this model the limits A — > a and — e — > — oo are 
determined while keeping the ratio of the second virial coefficient of the AHS model, -B AHS , 
to the second virial coefficient of the hard sphere system (HS), B^ s = 47rcr 3 /3, constant, 
that is 

oAHS 

%- = 4r AHS - 1 = const. (4) 

Equation (j4]) defines the stickiness parameter r AHS , which describes to which extend the AHS 
particles tend to glue together. Hence r AHS acts as an effective dimensionless temperature 
in the Baxter model. 

The Baxter model has been extensively investigated both by theory and by simulation.— - — 
Within the fluid phase regime it exhibits two types of phase transitions, namely a percolation 
transition and a vapor-liquid (vl) condensation transition. The critical parameters of the 
latter have been determined with very high accuracy using GCMC simulations. Their values 
are p A ^ s = 0.508(10)/a 3 (the critical particle density) and r A ^ s = 0.1133(5) (the critical 
temperature). Based on these values, we can now estimate the locus of the vapor-liquid 
critical point in the present square-well system. The first assumption is that the packing 
fraction rj = 7r<j 3 p/6 associated to the critical density remains constant, when we replace a, 
which corresponds to the average distance of two bonded particles in AHS system, with the 
corresponding average distance d av in the SW system. In other words, we require^ 1 ^ 

n J3 _!_ n AHS_3 (z\ 

For the present model, see Eq. (CE]), we have d av = cr + (A — cr)/2. Inserting this into Eq. (J5]), 
one obtains 

Pc - 1 ~ {a + {X-a)l2f (6) 



To estimate the critical temperature, we apply the Noro-Frenkel law of corresponding states. 
It states that two systems which are similar in nature by having both a hard-core repulsion 
and an (ultra-) short-ranged attraction have equal reduced second virial coefficients, provided 
we reduce with B^ s close to the critical temperatures in each case.— This allows us to apply 
a mapping of the phase diagram of one of those systems to the other. Taking the AHS as the 
reference system and equating its reduced second virial coefficient with the reduced second 
virial coefficient of a short-ranged SW model yields 

rAHS , rSW 

^hs- = 4r AHS -1^ = 1- ((A/a) 3 - l)(e^ - 1). (7) 

In Eq. ([7]) we introduced the reduced temperature T* = ksT/e. Solving Eq. ([7]) with respect 
to T* we obtain 
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In 1 + 



4r AHS ((2-A/a)- 3 - 1) J \ ' ® 
Based on Eqs. ([6]) and ([8]), we can estimate the critical parameters of the present model by 
inserting the known values for p^ s and t^ s , respectively. This yields 

77 CjV i = 0.250(5), T* vI = 0.346(5) 

where r] c ^\ is the critical packing fraction of the vapor-liquid condensation transition. 

In addition to determining the critical point, we may use the mapping procedure described 
above to obtain an estimate of the full fluid-fluid coexistence curve (binodal), as reviewed 



in Ref. 
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To this end, we apply Eq. (J6]) and Eq. ([8]) to GCMC data for the binodal taken 
from Fig. 3 of Ref. '3_1. Moreover, from integral equation studies^ of a SW system with 
A = 1.03a, the fluid-solid coexistence line is known (see Fig. 9 of Ref. 32). Taking their 
data and applying the mapping procedure again, we obtain an estimate of the fluid branch 
of the fluid-solid coexistence line in our system. Finally, we also mapped the data of the 



percolation threshold given in Ref. 
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Since we are partly far away from the critical point during such mapping, some doubts are 
appropriate regarding the justification of that strategy. Nevertheless, the mapping results 
give a first idea of the phase diagram of the present system. A summary of the results of 
these mapping procedures is given in Fig. [TJ Included in this figure are our present GCMC 
results for the bulk percolation thresholds, which will be discussed in more detail in Sec. IIVI 
The resulting percolation line divides the phase diagram into two regions characterized by 
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Figure 1. Estimated phase diagram of the bulk system, where we applied the mapping formulas © 
and ([8]) (empty symbols). O: the estimated critical point at r/ CjV i = 0.250(5) and T* vl = 0.346(5), 
A : the mapped metastable fluid- fluid coexistence line, o: the mapped percolation line and V: the 
mapped fluid branch of the fluid-crystal coexistence line. Our present results for the percolation 
threshold are indicated by filled symbols. All dotted lines are guides to the eye. 

finite (low rf) and infinite (high rf) clusters, respectively. Note that in both cases the clusters 
are transient in character, since there are no permanent bonds in our system. 

According to Fig. (TJ the lowest temperature T* = 0.5, for which we investigated the 
percolation transition, lies near or even within the metastable region between fluid and solid 
states. Therefore we carefully examined the packing fraction and configurations generated 
during each GCMC run in order to guarantee that we average only in one phase. In partic- 
ular, we analyzed the structure in terms of the radial distribution function g(r). In Fig. [2] 
we show g{r) together with the cavity function j(r) = g( y r)e~ Usw ^ r ^ kBT of the bulk system 
at T* = 0.5 at three packing fractions below, at and above the critical packing fraction for 
percolation rj c , which will be defined in Sec. IIVI Since g{r) and j(r) show no indication of 
any long rang order, we may conclude that we are indeed in the fluid phase. 

III. COMPUTATIONAL DETAILS 

To investigate the percolation behavior, we carried out standard GCMC simulations. 33 
For the bulk system calculations have been performed at several reduced temperatures 
T* = 0.5, 1.0, 1.5, 2.0. For the confined system we considered one fixed temperature, 
T* = 0.5, and several wall distances ranging from L z = 1.5a to L z = 10.0a. The sys- 
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Figure 2. a) Radial distribution function g{r) and b) cavity function y(r) at T* = 0.5 for packing 
fractions below, at and above the percolation threshold [rj ~ 0.1 (□), rj = rj c ~ 0.23 (O) and 77 ~ 
0.3 (A)]. To improve the visibility, the graphs for rj = r/ c {rj > rj c ) are shifted upwards by 1.0 (2.0) 
in a) and by 0.3 (0.6) in b). Dotted lines are guide to the eye. 

terns are equilibrated by performing at least 10 5 GCMC cycles, where each cycle consist on 
(N) single particle displacement, insertion or deletion attempts, with (JV) being the average 
number of particles in the system. Ensemble averages are obtained by analyzing at least 
10 5 equilibrium configurations, separated by 80 cycles. With this choice of parameters the 
energy autocorrelation function between two averaging events is close to zero in all cases. 

To determine the percolation threshold we performed at least four simulation runs at 
different basis lengths L for each system, where we have L x = L y = L z = L for the bulk 
system and L x = L y = L > L z for the slit-pore system. In each series, we chose the smallest 
L such that (N) « 500 at rj = 0.1. The largest systems we investigated contained about to 
30000 particles. 

An important aspect of any simulation study targeting percolation properties concerns 
the definition of clusters. Here we use a simple configurational criterion. Since we are 
dealing with a square-well model, two particles at positions T*j and rj are unambiguously 
defined as connected if their distance is within the attractive well of the pair potential, i. e., 
Tij < A. The resulting configurational clusters may be considered as "percolated" if they are 
connected to their own periodic image.— >^ In the 3D (bulk) system, we allow for percolation 
in any of the three spatial directions. The slit-pore system, on the other hand, is periodic 
only in the x- and the y-direction. Hence percolation can only occur in these two dimensions 
(2D). The whole system is considered as percolated if at least one percolating cluster exists. 
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To identify clusters and check for percolated configurations in our continuous model system, 
we implemented an algorithm as described in Ref. [15. 



IV. ANALYSIS OF THE PERCOLATION TRANSITION 
A. Determination of the percolation threshold 

The purpose of this section is to outline our approach to determine the percolation thresh- 
old. Indeed, the definition of percolation itself is not unique in the literature. It depends 
on the system to be investigated, the definition of a percolated state and the definition of 
bonds.- 1 ^ 1 ^ 5 - This may lead to incompatible results for a single system. On the other hand, 
the definition of the percolation threshold is much clearer, at least for lattice systems. There 
one typically considers the probability p of finding a given lattice site occupied. The thresh- 
old p c is then defined by the following criterion, that "for p above p c one percolating network 
exists; for p below p c no percolating network exists."— 

We now have to clarify which quantity in our continuous system should play the role 
of the control variable p. Since we perform GCMC simulations, a natural control variable 
would be the activity (or the chemical potential, respectively). Here we prefer to use the 
(average) packing fraction 77 due to the obvious analogy between 77 and the occupancy of 
lattice sites. 15 

To determine the percolation threshold r) c , we examine the probability H(rj, L) of finding 
a spanning cluster in a given system (3D or slit-pore) with a characteristic length L at 
packing fraction 77. This probability, which depends on the system size L and the packing 
fraction 77, can be calculated via an ensemble average of a test function 7r({rk}, L) 



where ir({rk}, L) — 1 if the configuration defined by {r^} has one spanning cluster and if 
all clusters are finite. It is commonly assumed (not rigorously proven) that 11(77, L) fulfills a 
single- variable scaling law^^ of the form 



The ansatz ( ITU]) involves the percolation threshold 7] c , the critical exponent v to be defined 
later [see Eq. (123]) ] and the master curve n(x), which is assumed to be unique for each set 



n(r],L) = (7r({r k },L)), 



(9) 



Tl(r ]j L) = fl(x) = fl((r ] -r ]c )L 1 ^). 



(10) 
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of parameters. To obtain an accurate estimate of i] c we employ the following strategy (for a 
more detailed description see Ref. [l). 

We first note that the probability H(r], L) is expected to have a sigmoidal shape with 
11(0, L) = and U(r] — > r] cp , L) = 1, where i] cp = /3 ~ 0.74 is the packing fraction at 
close-packing (cp). Given this shape, the first derivative of II(?7,L) with respect to i] is a 
function with a peak at the packing fraction characterizing the steepest ascent of U(r],L). 
This feature, combined with the fact that dU/drj is normalized to one, allows us to define 
an average density r^ av according to 

pep 

VM = j o d VV —. (11) 

Furthermore, we can define the standard deviation A av via 

rvcp f)r\ 
AavW = J o dv^-^) 2 —. (12) 

Using Eq. (TTU]) to substitute II for IT and writing 77 = x L~ l l u + r] c , we obtain from Eq. (ITT]) 
the linear dependence 

Vav = AL^ + Vc (13) 
where A is a constant. Applying similar arguments to Eq. (fT2|) results in 

A av oc L- l '». (14) 

Combining Eqs. (TLB"]) and (TT4"]) finally yields the expression 

77av - Vc OC A av . (15) 

Thus, Eq. (fl~5|) predicts a linear relationship between 7] av and A av , with r] c playing the role of 
an y-intercept. We have used this relationship to obtain the percolation threshold. In order 
to actually calculate ?7 av and A av , we have fitted the numerical data for Il(?7, L) according 



to the expression (see Ref. Il5l ) 
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/(x;{a i }) = l + tanh^o i x i J (16) 

involving the six fitting parameters {a{\. Some representative numerical results for ^ av 
as function of A av are given in Fig. [3^) and c). In Fig. [3^) we consider bulk systems at 
various reduced temperatures. Clearly, the data points for r/ av form nearly perfect straight 
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Figure 3. Upper panel: Determination of the percolation threshold and the exponent v for bulk 
systems of temperatures T* = 0.5 (•), 1.0 (A), 1.5 (T) and 2.0 (♦). The plots show a) the average 
packing fraction ry av as function of the standard deviation A av and b) the (logarithm of) the standard 
deviation as function of system size. The dashed lines are obtained from least square fit. From 
a) the percolation threshold follows as the intercept on the y-axis [see Eq. (|15p j. From b), the 
exponent v follows as the slope of the lines [see Eq. (|14j) ]. Lower panel: Analog to the upper panel, 
but for slit-spore systems at temperature T* = 0.5 and L z = 3. Oct (•), 3.5o" (A) and 4.0a (T). 



lines, consistent with Eq. (|T5l) . which allows us to extract the percolation threshold density 
T] c easily. The averaged gradient A of these nearly parallel lines is essentially zero, A = 
—0.002 ± 0.005. The results for rj c at the four temperatures investigated are plotted in 
Fig. [TJ Figure |3fc) shows the same type of data for the slit-pore systems with various wall 
separations. The results reflect a very good accuracy of the approach again. Here, the 
averaged gradient is close to —1/2 (i. e., A = —0.505 ± 0.005). The resulting percolation 
thresholds for the confined systems are further discussed in Sec. IIVBI 

In a similar fashion, we can use Eq. f fl^|) to obtain the critical exponent v [which is involved 
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in the scaling law (ITU]) ]. To this end we plot in Fig. [3b) and d) the standard deviation A av 
as a function of the system size in a double-logarithmic representation. Again, we find that 
the relationship predicted by Eq. (fl4"]) is essentially perfectly fulfilled for both, bulk and 
slit-pore systems. Finally, having determined r] c and v, we can test the underlying scaling 
assumption for Il(?7, L) given in Eq. (1 10 1) . We have done this test for all systems investigated. 
Results for bulk and slit-pore systems at T* = 0.5 are shown in Fig. HJ For both cases the 
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Figure 4. Upper panel: Percolation probability H(rj, L) of the bulk system at T* = 0.5. a) Simu- 
lation data for three system sizes L = 8. 06a (□), 16.12<7 (O) an d 29.7a (a). The dashed lines are 
obtained via least square fits of Eq. (TLB"]) . The bold solid line denotes the assumed limiting shape 
for L — > oo. The dotted line marks Il(r/ C , L). b) Shifted and scaled simulation data (x) for all nine 
investigated system sizes, ranging between L = 8.06<r and L = 29.7a. The bold solid line denote 
the assumed shape of the master curve 11(2;) as obtained from a least square fit of Eq. (|16p . Lower 
panel: Same as the upper panel, but for a slit-pore system at T* = 0.5, L z = 3.0a and three system 
sizes L = 13.21cr (□), 29.54a (O) and 93.42cr (a). 

data obtained for different system sizes indeed falls on one master curve, consistent with the 
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theory. Numerical values for the exponent v are given in the subsequent Sec. IIV B[ 



B. Further critical exponents 

In addition to the exponent v mentioned in Sec. IIV Al we have determined three more 
critical exponents. These are the two Fisher exponents (FE) tfe and ctfe ; and the fractal 
dimensionality D, which are related to the cluster size distribution, n(s;r]), and the distri- 
bution of radii of gyration, R(s; 7/). 4,37 Thus, the exponents tfe, Cfe, and D have geometric 
(rather than thermodynamic) character. Nevertheless, one can show- that these exponents 
are linked to thermodynamic quantities such as the pressure via scaling relations. This link 
can be made in the framework of the fisher droplet model 3 ^" 3 ^ of an ideal cluster gas. Here 
we focus on relationships between the "geometric" exponents. The cluster size distribution 
is defined as the average number of finite clusters of size s at packing fraction rj. Close to 
the percolation threshold, this quantity displays power-law behavior, i. e., 

n(s; Vc) ~ s~ TFE (17) 

as s — > oo. Equation fl!7j) defines the Fisher exponent tee- The second Fisher exponent, 
apE; appears when we consider the k-th moment m& of n(s;r)) 

oo 

m k (r]) = ^2s n n(s;7]). (18) 

s=0 

Using Eq. (TTT|) . we find that m k fulfills the power law 

m k (r)) oc \rj - r/ c |- (fe - T+1)/fJFE . (19) 

We now consider the distribution of radii of gyration R(s; 77). For a fixed size s, the radius of 
gyration of a single cluster is defined as the averaged mean squared distance of the particles 
relative to the cluster's center of mass (CM) position t*cm — Y2i=i r i/ s - This implies^ 




Close to the percolation threshold, the critical behavior of this distribution is given through 
the power law 

#McW /D (21) 
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as s — > oo. We introduced the fractal dimensionality D as the critical exponent of the 
gyration radius. The distribution R(s; rj) also allows us to define the correlation length 
characterizing the system, 



|2 Em«AH. (22) 

The correlation length diverges at the percolation threshold as^ 1 ^ 

£ oc |tj - Tjcr, (23) 

which defines the exponent i/ already mentioned in Sec. IIV A[ Finally, using the definition 
(122]) and the power- law relations f fT9|) . (12T|) and (123]) one can find the scaling relation 

£> = (24) 
a FE u 

which we use to determine the fisher exponent <7fe from the knowledge of v and D. The 
power laws given in Eqs. [T7] and [21] hold only in the limit s — > oo. In fact, for the cluster 
size distribution the exponent Tfe determines only the leading term in a series expansion of 
the form 

n(s;f) c )~r FE (l + Ar !1 + ...). (25) 

In Eq. (125]) the critical exponent Q is introduced as a first order correction to the scaling 
ansatz.— 

To treat the system of finite size, we performed simulations for each analyzed slit-pore 
system, using the appropriate values of the critical activity to obtain r] c in each case. We 
prepared the raw data by "logarithmic binning" of the distributions n(s; r] c ) and R(s; r] c ). To 
this end, we average the distributions over intervals I n = [feo 6 "^" 1 , bo ^T 1 ], where the width 
b n of each /„ grows as b n = bob n and hence produces equally spaced nodes on a logarithmic 
axis. — This introduces some arbitrariness via the choice of b. For practical reasons we have 
to choose b such that the noise-level of the averaged distributions is sufficiently small. On 
the other hand, we must be careful not to obscure the main tendencies. We then plot the 
averaged data in a double-logarithmic representation and approximate the local slope by 
difference quotients. For sufficiently large s we expect to observe a clear plateau. Indeed, 
from the definition of the critical exponents rpg and D in Eqs. (T2~T]) and (125]) it follows that 

lim n'(s; r] c ) = -r FE (26) 
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and analogously 

\imR'(s;r] c ) = ±-. (27) 

s— >oo J_J 

Here, f'(s) = din f(s)/ dins (where / = n or R). In Fig. [5] we show exemplary data for the 
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Figure 5. Local slope n'(s;rj c ) of the double-logarithmic representation of the cluster-size distribu- 
tion n(s; rj c ) (inset) for the bulk system at reduced temperature T* = 0.5 and four system sizes 
L = 9.35a (A), 18.71a (•), 37.42a (o) and 55.13a (■). The logarithmic-binning (See Sec. ITVBjl is 
done using b = 2 for the local slope and b = 1.3 for the inset. The dashed line denotes the value of 
the fisher exponent tfe for 3D random percolation. 

functions n'(s; rj c ) and n(s; r] c ) (inset) for the bulk system at T* = 0.5 for four system sizes. 
For the largest system, L = 55.13a, we can indeed clearly identify a plateau. The deviations 
for large s are due to the limitations in the system size. For clusters with a characteristic 
length approaching the system length L, the probability of being not percolated has to drop 
rapidly, since percolated clusters do not contribute to the distributions. 

In Fig. |6] we show, again for the bulk system, the two functions n'(s;r] c ) and R\s;r] c ) 
at four different temperatures (and the largest L considered). All curves exhibit a plateau, 
which is reached the faster (in terms of the cluster size s) the higher T*. The important 
point, however, is that the plateau heights are essentially independent of T*. This means 
that the exponents tfe and D [see Eqs. fl26l) and ( 1271) ] extracted from our simulations may 
be considered as "universal", that is, not influenced by the temperature. Furthermore, the 
observation that the temperature does not affect the value of the critical exponents but the 
speed of convergence may offer an explanation for the observed temperature dependence of 
the critical exponents for more complex interaction models.— An overview of our results for 
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Figure 6. Logarithmic slopes n'(s; r/ c ) a) and R'(s; rj c ) b) for the bulk system at the four investigated 
temperatures T* = 0.5 (o). 1.0(»), 1.5(A) and 2.0(A) directly at the corresponding percolation 
thresholds. The results are for the largest system size investigated for the bulk system (L = 55.13cr). 
The bold dashed lines denote the values for the critical exponent tfe (part a)) and D (part b)) for 
3D random percolation. 

the bulk system is given in Tab. [H It is visible that the critical exponents tfe and D are 



Table I. Critical parameters for the bulk system. 



Vc 



Il c V D TFE <7FE 



0.5 0.23165(5) 0.47 0.86 2.48 2.19 0.47 

1.0 0.32960(5) 0.44 0.87 2.51 2.20 0.46 

1.5 0.36337(5) 0.46 0.86 2.52 2.21 0.46 

2.0 0.3800(3) 0.48 0.88 2.51 2.21 0.45 

3D r. p.* - - 0.88 2.53 2.18 0.45 



a Literature values for the critical exponents relevant in 3D random percolation (r. p.)ild2 

close to those obtained for 3D random percolation (see horizontal lines in Fig. [BJ. 

We now turn to the slit-pore systems. In Fig. [7] we present n'(s;r) c ) and R'(s;r] c ) as 
functions of ln(s) for six values of the wall separation L z . The convergence of the data (as 
characterized by the appearance of a plateau) strongly depends on L z . Specifically, the 
larger L z , the larger cluster sizes are needed before n' and R' reach a plateau. Moreover, in 
some cases, there is no clear plateau at all (see e.g. data for n' at L z = 10a). For the other 
cases, the observed plateau heights are close to those predicted by 2D random percolation 
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Figure 7. Same as Fig. [6l but for slit-pore systems at temperature T* = 0.5. To improve the 
visibility, only the results for L z = 1.5a (A), 3.0<t(a), 4.0a (y), 5.0a (T), 7.5a (0), and 10.0a (♦) 
are shown. The results are obtained for the largest system size we investigated in each case. The 
bold dashed lines denote the values for the critical exponent r (part a)) and D (part b)) for 2D 
random percolation. 

theory. The latter are summarized in Table [Til An overview of our numerical data for the 
Table II. Critical exponents for 2D random percolation. 





l/ 2D 


D 2D 


r 2D 
FE 


rr 2D 
°FE 


exact 


4 
3 


91 
18 


187 
91 


36 
91 



numerical 1.33 1.90 2.05 0.40 



exponents tfe, D, u, and cp^ as functions of the wall separation is given in Fig. [HJ The 
graphs also indicate the corresponding literature values for 2D and 3D random percolation. 
Within the entire range investigated (L z < 10a) the fractal dimensionality is very close to 
the 2D value, and a similar consistency is found for the exponent of the correlation length, 
v, and the exponent aFE- Compared to that, the data for tfe are more scattered, consistent 
with what is seen in Fig. [7^). 

As far as the slit-pore systems are considered, there are (at least) two sources of errors 
which come into play. The larger the wall separation L z , the larger cluster sizes s are 
needed to observe a clear plateau. Hence, our system size may still be chosen too small. 
Furthermore, we can observe real 2D behavior only if we know the percolation thresholds 
T) c with sufficient accuracy. We expect the error tolerance of r] c to decrease with increasing 
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Figure 8. Estimated critical exponents tfe ( a ); D (b), v (c) and <7fe (d) for the slit-pore systems 
at T* = 0.5 as functions of the wall distances L z . The exponent ctfe is deduced from Eq. (|24p . The 
lines denote the values for 3D (dashed) and 2D (dotted) random percolation. 

wall distance L^M^ and in this spirit, we only measure effective exponents. At L z w 2a 
we observe, especially for the fisher exponent tfe, strong deviations from the 2D value of 
random percolation. We suppose this irregularity to occur due to the vicinity a freezing 
transition, which position in the phase diagram is known to depend on the wall separation^ 
(see Fig. [9] a). 



C. Confinement-induced shift of the percolation threshold 

So far we have focused on the critical exponents characterizing the percolation transition. 
We now consider, for the confined systems, the actual location of this transition in the phase 
diagram. To this end, we depict in Fig. 0a) the critical packing fraction r\ c and in Fig. 0b) the 
corresponding chemical potential of the percolation threshold for various wall separations. 
For comparison, we also show the freezing line of a pure hard sphere system confined between 
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Figure 9. Relation between the percolation threshold and the wall separation for the slit-pore 
systems at T* = 0.5. Specifically, part a) shows the packing fraction r] c (L z ), while b) shows the 
difference of the chemical potential to the corresponding bulk value /3A/i c = f3(p, c (L z ) — /x c (oo)). 
The dotted line is a guide to the eye. The dashed line denotes the value of the bulk threshold at 
rj c (oo) [see Tab. HI] . The solid line marks the freezing line of a pure hard sphere system which is 
confined between two hard walls (data taken from Fig. 4 of Ref. |45| ). 



two parallel, hard walls (solid line) in Fig. [9] a). The data was taken from Fig. 4 of Ref. |45 
The pure hard-sphere system corresponds to our model at large temperatures, i. e., T* — > oo. 

Inspecting Fig. [9]a) we find that with decreasing L z the percolation threshold i] c increases 
until L z « 2er. At even smaller wall separations the relation between rj c and L z becomes 
non-monotonic. 

In the following we focus on the behavior observed for L z > 2a, i. e., the shift of the 
percolation threshold with respect to the bulk values. In particular, we give an argument 
which has successfully been used to explain the shifts of the gas-liquid critical point in a slit- 
pore confinement.— 1 ^ The most remarkable point of this argument is that it does not rely on 
any detailed information of the underlying system of interest, but only on the dimensionality 
and the interaction range. In Sec. IIV Al we have derived the percolation thresholds for fixed 
L z by performing a finite-size analysis with respect to the lateral dimensions L x = L y = L 
[see Eq. (fTTj) ]. i. e., we have let L — > oo. This way we obtaine the percolation threshold r] c , 
which still depends on the remaining length L z . We now see that, for large L z , the behavior 
upon further increasing L z towards infinity is the same as that seen for L — > oo for the bulk 
system.— In other words, we assume that a relation similar to Eq. (TTTj) should also hold for 
rj c (L z ). Moreover, we assume that, for sufficiently large L z , the corresponding exponent v 
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should be that of the 3D system, v m 0.88. These considerations lead to 

A V * C (L Z ) = r]*(L z ) - V * c (oc) cc L- 1 '^. (28) 

In Eq. (128|) . r]*(L z ) denotes the packing fraction of a bulk system which is in chemical 
equilibrium with a slit-pore system of width L z (and at average packing fraction r] c ). Indeed, 
using this notation, rj* (oo) is identical to the percolation threshold of the bulk system rj^ D 
at T* = 0.5. 

Since we consider our system to be above the gas-liquid critical temperature the density 
is a continuous function of the chemical potential. Hence we may write, for small deviations 
Afi c {L z ) from the bulk value /i c (oo) (i. e., sufficiently large L z ), 

A/i c (L 2 ) = fi B (L g ) - fi e (oo) oc L~ 1/um . (29) 

In Fig. [10] we compare our numerical results for Ar]*(L z ) and Afi c (L z ) with the predictions 
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Figure 10. Scaling of the confinement-induced shifts of the percolation threshold at T* = 0.5. The 
shifts are given in terms of a) the corresponding densities of a bulk system in equilibrium with the 
slit-pore Arj*, and b) the change of the chemical potential /3A/i c . The gradient of the straight lines 
is set to —1/vzd as required by Eqs. (I28p and (|29p . 

of Eqs. (I28|) and f )29|) . by plotting a straight line (dashed) with a fixed slop — 1/vzn together 
with the data. From Fig. [10] a) we find that Eq. ([28]) is at least qualitatively fulfilled even 
for the smallest wall distances investigated. Moreover, the shift of the chemical potential 
shown in Fig. [10] b) is sufficiently well described by Eq. (|29|) for L z > 4a. For not too small 
L z we can thus conclude that the confinement- induced shift of the percolation threshold r) c 
can successfully be described by Eqs. ([28]) and ([29]) . Hence we have shown that the same 
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scaling arguments, which are usually applied to the shift of the gas-liquid critical point in a 
slit-pore confinement,— 1 ^ also apply to the shift of the percolation threshold. 



V. CONCLUSION 

Based on grand-canonical MC simulations and a finite-size-scaling analysis we have stud- 
ied the percolation behavior of a model colloidal suspension with an ultra-short range of at- 
traction in three dimensions and in slit-pore geometries. Our aim in this study was twofold: 
On the one hand, we wanted to provide "fresh" data on the somewhat controversial issue of 
the universality of critical exponents of percolation in 3D continuous systems. On the other 
hand, we aimed to explore the impact of spatial confinement on the percolation transition; 
a question which was so far essentially unexplored. The model we have studied consists of 
hard spheres with a square-well attraction extending over (only) four percent of a particle 
diameter. As a result, the bulk systems lacks a stable liquid phase (i. e., the vapor-liquid- 
transition predicted by a mapping procedure from the Baxter model is metastable), and the 
only true transition within the fluid phase is the percolation transition. Using a finite size 
analysis we have first determined the percolation threshold, r] c , of the bulk system at four 
temperatures above the critical temperature. Consistent with an earlier simulation study of 
a similar square- well fluid (yet with a significantly larger attraction range) ,— we found that 
T) c increases with increasing reduced temperature and slowly converges to the limiting value 
determined by a pure hard sphere system (i. e., the limit T* — > oo of our model). 

Regarding the critical exponents of the 3D percolation transition, however, our present 



results differ from the earlier ones in Ref. [14|. In that study, the exponent v of the correlation 
length was found to depend on the temperature, specifically v decreased with T*. This led 
the authors of Ref. 14 to conclude that the exponent is "non-universal". Contrary to these 
results, our present data reveal no systematic dependence of v on the temperature; instead, 
the numerical values of v fluctuate around the corresponding value predicted by 3D random 
percolation theory. The question remains whether this contradiction could be due to the 
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different range of attraction considered in our study and in Ref. 

In addition to we also studied the exponents tfe and D governing the cluster size 
distribution n(s; r] c ) and the radii of gyration, respectively. For sufficiently large system 
sizes, we find these exponents to be independent on the temperature and their magnitude 
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to be very close to those predicted by 3D random percolation theory again. 

To get an insight of the impact of spatial confinement, we performed GCMC simulations 
at several wall separations in the range 1.5a < L z < 10a at fixed reduced temperature 
T* = 0.5. All of these confined systems displayed a clear percolation transition, with the 
percolation value of the chemical potential (and also the resulting r/ c ) being shifted relative 
to its bulk value (at T* = 0.5). 

Specifically, upon decreasing L z from infinity (bulk limit), we find a monotonic increase of 
r) c until the wall separation comes close to the 2D limit {L z j o > 2) . Moreover, in this range of 
wall separations the shift of r] c can be mathematically described by simple scaling arguments 
similar to those applied in the framework of cross-over scaling for vapor-liquid critical points 
in slit pores.— ^ In addition to the percolation threshold, we have also analyzed the critical 
exponents for the slit-pore systems (see Fig. [8]). The exponent v is always close to the value 
predicted by 2D random percolation theory and shows no systematic dependence on L z . 
Regarding the other exponents, we have found that convergent results (extractable from a 
plateau in the derivatives of the relevant distribution functions) are generally the harder to 
obtain the larger the wall separation is. Nevertheless, in the cases where we could extract 
reliable data, the exponents are comparable to those of 2D random percolation. 

We thus conclude that the percolation transition of colloidal systems with ultra-short 
ranged interactions falls into the universality class of random percolation, the only important 
factor being that the spatial dimension along which the system is finite (3D or 2D). 

From an experimental and applicational point of view one is often interested in the 
consequences of percolation, such as a strong increase of the macroscopic conductivity, rather 
than in the transition itself. In this context, we expect our results for the confined systems 
to be particularly relevant, since real experimental situations often involve confinement. We 
note that the model considered here is appropriate for a broad range of colloidal suspensions 
involving depletion "agents" which cause the short-range attraction. From a more conceptual 
point of view, it is interesting that colloidal systems with such "sticky" interactions — at 
least in 3D — can not only form percolated static structures, but an also display non-trivial 
dynamical behavior such as gelation and glass formation.— In particular, for spatially 
confined systems the interplay between gelation and percolation is essentially unexplored. 
These issues will be subject of future studies. 
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